doseresNMA antidep/check knots .R

In htx-r/doseresNMA:

# this file is to run the model where knots are generated from the dose range 0:max to compare with min:max
# libraries
library(R2jags)

# functions
source('data.prep.R') # prepare data
source('net.structure0knot.R') # jags model format
source('drnma.model.standardRE.R') # standard RE jags model

# data
antidep <- data.prep()

#-- Model1 - standard RE delta's --
# data
jagsdata_std_0knot <- net.structure(data=antidep)
# run model
drnma_RE_0knot <- jags.parallel(data = jagsdata_std_0knot,inits=NULL,
                          parameters.to.save = c('beta1','beta2','g','tau.d','p.drug','resdev','rhat','totresdev','Z','r','n'),
                          model.file = drnma.standardRE,
                          n.chains=3,n.iter = 10000,n.burnin = 4000,n.thin = 1,
                          DIC=F)
save(drnma_RE_0knot,file='drnma_RE_0knot')

# plot
load('drnma_year')
source('drnma.plot.2absolute.R')
ggsave('fig2_0knot.pdf',
       absolutePredplot2(x=drnma_RE_0knot),
       width =8.27,
       height =11.69)

# create table with both knots: start from 0 and min.dose
# 1. take the drug-dose and find its RCS transformation at
dose_to_rcs_min <- function(dose.per.drug,knot_probs=c(0.25,0.50,0.75)){
  require('rms')
  max.dose <- max(dose.per.drug)
  min.dose <- min(dose.per.drug)
  knots <- quantile(min.dose:max.dose,probs = knot_probs)
  rcs.dose.per.drug <-  rcspline.eval(dose.per.drug,knots = knots,inclx = TRUE)
  return(knots)
}


dose_to_rcs_0 <- function(dose.per.drug,knot_probs=c(0.25,0.50,0.75)){
  require('rms')
  max.dose <- max(dose.per.drug)
  min.dose <- min(dose.per.drug)
  knots <- quantile(0:max.dose,probs = knot_probs)
  rcs.dose.per.drug <-  rcspline.eval(dose.per.drug,knots = knots,inclx = TRUE)
  return(knots)
}

# min.doses
knots_min <- antidep[antidep$drug!='placebo',] %>%
  group_by(drug) %>%
  group_map(~dose_to_rcs_min(.x$dose)) # the list order as the order of levels(data_no_placebo)
tbl_min <- do.call(rbind,knots_min)
# 0 doses

knots_0 <- antidep[antidep$drug!='placebo',] %>%
  group_by(drug) %>%
  group_map(~dose_to_rcs_0(.x$dose)) # the list order as the order of levels(data_no_placebo)
tbl_0 <-do.call(rbind,knots_0)

# save 
tbl_min_0 <- cbind(tbl_min,tbl_0)
write.table(tbl_min_0,'knots_min_0')
cbind(tbl_min_0[,4])